#define SHADE

const float OVERFLOW = 1.5e38;

const int N=100, M=200;

const vec4 c=vec4(-0.08,0.0,-0.8,-0.03);

const vec3 light0=vec3(1.0, 1.0, 1.0);
const vec3 light1=vec3(-.5, 0.866025, 1.0);
const vec3 light2=vec3(-.5, -0.866025, 1.0);
const vec3 light3=vec3(0, 0, 1.0);

const vec3 lightColor0=vec3(.25, .0, .0);
const vec3 lightColor1=vec3(.0, .25, .0);
const vec3 lightColor2=vec3(.0, .0, .25);
const vec3 lightColor3=vec3(.8, .8, .8);

uniform float w;
uniform float epsilon;
uniform mat4 transform;

vec4 qmult(vec4 q1, vec4 q2) {
	return vec4(
		q1.x*q2.x - dot( q1.yzw, q2.yzw ),
		q1.x*q2.yzw + q2.x*q1.yzw + cross( q1.yzw, q2.yzw )
	);
}

vec4 q2(vec4 q) {
	return vec4(
		q.x*q.x - dot( q.yzw, q.yzw ),
		2.0*q.x*q.yzw
	);
}

float d(vec4 x, float epsilon) {
	vec4 dx;
	float lx2;
	int i;
	
	dx=vec4(1.0, .0, .0, .0);
	
	for(i=0; i<N; i++) {
		lx2=dot(x, x);
		if(lx2 > 3.0) 
			break;
		
		dx=2.0 * qmult(x, dx);
		x=q2(x) + c;
	}
	
	if(lx2 < 3.0) {
		return 0.0;
	} else {
		float ldx2=dot(dx, dx);
		
		if(ldx2 > OVERFLOW) {
			return .0;
		} else {
			return .25 * sqrt(lx2 / ldx2) * log(lx2);
		}
	} 
}

float lshade(vec3 n, vec3 l, vec4 x) {
	return max(.0, dot(n, normalize(l - vec3(x))));
}

void color(vec4 x, float bounding) {
	float epsilon=bounding / 2.0;

	vec4 xp=vec4(epsilon, .0, .0, .0);
	vec4 yp=vec4(.0, epsilon, .0, .0);
	vec4 zp=vec4(.0, .0, epsilon, .0);
	
	vec3 n=vec3(
		d(x + xp, epsilon) - d(x - xp, epsilon),
		d(x + yp, epsilon) - d(x - yp, epsilon),
		d(x + zp, epsilon) - d(x - zp, epsilon)
	);
	
	#ifdef SHADE
	n = normalize(n);
	
	vec3 lr, color;
	
	color 
		= vec3(.02)
		+ (lightColor0 * lshade(n, light0, x))
		+ (lightColor1 * lshade(n, light1, x))
		+ (lightColor2 * lshade(n, light2, x))
		+ (lightColor3 * lshade(n, light3, x))
		;
	
	gl_FragColor=vec4(color, 1.0);
	#else
	gl_FragColor=vec4(.9);
	#endif
}

void main()
{
	vec4 s, x;
	float delta, z;
	int i;
	
	s=gl_TexCoord[0];
	
	gl_FragColor=vec4(.0, .0, .0, 1.0);
	
	z=2.0;
	i=0;
	while(length(z) < 4.0) {
		x=vec4(s.x, s.y, z, w);
		delta=d(transform * x, epsilon);
		
		if(delta == 0.0) {
			color(x, epsilon);
			break;
		} else {
			if(delta <= epsilon || i > M) {
				color(x, delta);
				break;
			} else {
				z -= delta;
			}
		}
		
		i++;
	}
}
